Figures

1 Import library

Code

# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)
library(gridExtra)
library(png)

# table
library(gt)
library(gtable)

# map
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle and timezone
library(suncalc)
library(lutz)

# data wrangling
library(tidyr)
library(dplyr)
library(data.table)
library(purrr)
library(forcats)
library(lubridate)

2 Setting up custom function

2.1 matlab_to_posix

Code

# function to convert matlab date into R posix
matlab_to_posix <- function(x, timez = "UTC") {
  days <- x - 719529 # 719529 = days from 1-1-0000 to 1-1-1970
  secs <- days * 86400 # 86400 seconds in a day
  return(as.POSIXct(
    strftime(
      as.POSIXct(secs,
        origin = "1970-1-1",
        tz = "UTC"
      ),
      format = "%Y-%m-%d %H:%M:%S",
      tz = "UTC",
      usetz = FALSE
    ),
    tz = timez
  ))
}

2.2 windrose

Code

windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet_wrap = F,
           facet_grid = F,
           hour_sunset = 20,
           hour_sunrise = 7) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, hour_sunrise, hour_sunset), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin labels
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet_wrap == "origin") {
      data2$origin <- unique(data_to_plot$origin)
    } else if (facet_wrap == "trip") {
      data2$trip <- unique(data_to_plot$trip)
    }

    # add trip and origin if any
    if (facet_grid) {
      data2$origin <- unique(data_to_plot$origin)
      data2$trip <- unique(data_to_plot$trip)
    }

    # set up max value
    maxvalue <- 15

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "black",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      # somehow was needed
      theme(panel.grid.major = element_blank()) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet_wrap == "origin") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    } else if (facet_wrap == "trip") {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(trip ~ .,
          strip.position = "right"
        )
    }

    # if facet
    if (facet_grid) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_grid2(trip ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code

# import the data
# data_dive <- readRDS("../export/data_dive_rev_3_inter.rds")
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartLoc == "ANNU")

# remove id == `2015006` (no benthic dives)
data_dive <- data_dive %>%
  filter(!id %in% c("2015006"))

Code

# get predator distribution *.kmz (extraction using google earth)
# for calculation
Shark_area_coord <- read_sf("../export/shark.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  dplyr::select(lon = X, lat = Y)
Orca_area_coord <- read_sf("../export/orca.kmz") %>%
  smooth(., method = "ksmooth", smoothness = 0.5) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  dplyr::select(lon = X, lat = Y)

# for plots
Overlap_area_coord_to_display <- st_intersection(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  st_coordinates() %>%
  as_tibble() %>%
  dplyr::select(lon = X, lat = Y)
Shark_area_coord_to_display <- st_difference(
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  dplyr::select(lon = X, lat = Y)
Orca_area_coord_to_display <- st_difference(
  read_sf("../export/orca.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5),
  read_sf("../export/shark.kmz") %>%
    smooth(., method = "ksmooth", smoothness = 0.5)
) %>%
  # st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%
  st_coordinates() %>%
  as_tibble() %>%
  dplyr::select(lon = X, lat = Y)

Code

# # change few things
# data_dive <- data_dive %>%
#   mutate(
#     # change divetype
#     DiveType = if_else(DiveType == 3 &
#                          !is.na(Lat) &
#                          (abs(bathy) - Maxdepth) > 100, 0, DiveType),
#     # create DiveTypeName
#     DiveTypeName = case_when(
#       DiveType == 0 ~ "Transit",
#       DiveType == 1 ~ "Forage",
#       DiveType == 2 ~ "Drift",
#       DiveType == 3 ~ "Benthic"
#     ),
#     # add Post Molting vs Post Breeding trip
#     trip = if_else(trip == "PB", "Post Breeding", "Post Molting"),
#     # add date
#     date = as.Date(DateTime),
#     # add inside/outside each predators area
#     shark_area = point.in.polygon(
#       point.x = Lon,
#       point.y = Lat,
#       pol.x = Shark_area_coord$lon,
#       pol.y = Shark_area_coord$lat
#     ),
#     orca_area = point.in.polygon(
#       point.x = Lon,
#       point.y = Lat,
#       pol.x = Orca_area_coord$lon,
#       pol.y = Orca_area_coord$lat
#     ),
#     overlap_area = as.numeric(if_else(shark_area == 1 &
#                                         orca_area == 1, T, F)),
#     total_area = as.numeric(if_else(shark_area == 0 &
#                                      orca_area == 0, F, T))
#   )

4 Data Visualisation

4.1 Table 1

We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.

Code

compute_table_1_a <- function(x) {
  # if the sum of the column is below 0.1
  if (sum(x) < (0.1 * 2)) {
    # it's the proportion of dives that are benthic
    # so we weighted average these by taking into account all dives
    y <- weighted.mean(x, c(
      data_dive %>% filter(trip == "Post Breeding") %>% nrow(),
      data_dive %>% filter(trip == "Post Molting") %>% nrow()
    ))
    y <- paste0(round(y * 100, 1), "%")
    return(y)
    # otherwise if below 1
  } else if (sum(x) < (1 * 2)) {
    # it's the proportion of benthic dives within predator areas
    # so we weighted average these by taking into account only dives with location data
    y <- weighted.mean(x, c(
      data_dive %>% filter(trip == "Post Breeding" & !is.na(Lat)) %>% nrow(),
      data_dive %>% filter(trip == "Post Molting" & !is.na(Lat)) %>% nrow()
    ))
    y <- paste0(round(y * 100, 1), "%")
    return(y)
    # otherwise we just sum the column
  } else {
    y <- prettyNum(sum(x), big.mark = ",")
    return(y)
  }
}

Code

# build table
table_1_a <- data_dive %>%
  group_by(trip, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) /
      sum(n_benthic_dives_w_loc, na.rm = T)
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives per area",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark",
    n_benthic_dives_in_orca_proportion = "Orca",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_a")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  # color rows
  opt_row_striping() %>%
  # width
  tab_options(table.width = pct(130))

Code

# custom function to summarize each column differently
compute_table_1_b <- function(x) {
  # if the column is numeric
  if (all(is.numeric(x))) {
    # and the sum is below 100
    if (sum(x) < 100) {
      # it is time
      # so we "circular" weighted average them
      y <- psych::circadian.mean(
        c(
          rep(x[1], data_dive %>%
            filter(trip == "Post Breeding" & !is.na(Lat)) %>%
            pull(id) %>%
            n_distinct()),
          rep(x[2], data_dive %>%
            filter(trip == "Post Molting" & !is.na(Lat)) %>%
            pull(id) %>%
            n_distinct())
        ),
        na.rm = T
      )
      y <- round(y, 1)
    } else {
      # otherwise it is just number
      # so we "just" weighted average them
      y <- weighted.mean(x, c(
        data_dive %>%
          filter(trip == "Post Breeding") %>%
          pull(id) %>%
          n_distinct(),
        data_dive %>%
          filter(trip == "Post Molting") %>%
          pull(id) %>%
          n_distinct()
      ))
      y <- round(y, 1)
    }
  } else {
    y <- "\u00b1"
  }
  return(y)
}

Code

# build table
table_1_b <- data_dive %>%
  group_by(trip, id) %>%
  group_by(trip, id) %>%
  # create a column to say when start and end each trip by id
  mutate(
    is_arrival = if_else(DiveNumber == max(DiveNumber), T, F),
    is_departure = if_else(DiveNumber == min(DiveNumber), T, F)
  ) %>%
  # then summarise
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = hour(date_tz[is_departure == T]) +
      minute(date_tz[is_departure == T]) / 60 +
      second(date_tz[is_departure == T]) / (60 * 60),
    hour_day_arrival = hour(date_tz[is_arrival == T]) +
      minute(date_tz[is_arrival == T]) / 60 +
      second(date_tz[is_arrival == T]) / (60 * 60),
    .groups = "drop"
  ) %>%
  group_by(trip) %>%
  summarise(
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt(rowname_col = "trip") %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # rename columns
  cols_label(
    trip = "Trip",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # summary rows
  grand_summary_rows(
    columns = everything(),
    fns = list(label = "OVERALL", fn = "compute_table_1_b")
  ) %>%
  # color rows
  opt_row_striping()

Code

# display first table
table_1_a

Table 1:

Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

	# seals	All dives			Benthic dives		% Benthic dives per area
		Total	w/ loc	w/o loc	Total	Proportion	Shark	Orca	Total	Overlap
Post Breeding	244	1,090,092	972,017	118,075	33,982	3.1%	67.5%	83.0%	83.2%	67.2%
Post Molting	233	2,597,260	1,980,478	616,782	51,748	2.0%	70.4%	87.2%	87.4%	70.2%
OVERALL	477	3,687,352	2,952,495	734,857	85,730	2.3%	69.4%	85.8%	86%	69.2%

Code

# display second table
table_1_b

Table 2:

Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

	# dives			# benthic dives			Hour departure			Hour arrival
	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
Post Breeding	4467.6	±	964.6	139.3	±	230.7	20.6	±	1.7	9.3	±	2.8
Post Molting	11147.0	±	2093.1	222.1	±	384.3	19.3	±	1.7	13.7	±	1.9
OVERALL	7730.3	±	1515.8	179.7	±	305.7	20	±	1.7	11.2	±	2.4

Notes

• Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.

• Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals

Code

# get the hours of departure and arrival
hour_data <- data_dive %>%
  group_by(trip, id) %>%
  # create a column to say when start and end each trip by id
  mutate(
    is_arrival = if_else(DiveNumber == max(DiveNumber), T, F),
    is_departure = if_else(DiveNumber == min(DiveNumber), T, F)
  ) %>%
  # test
  summarise(
    hour_day_departure = hour(date_tz[is_departure == T]) +
      minute(date_tz[is_departure == T]) / 60 +
      second(date_tz[is_departure == T]) / (60 * 60),
    hour_day_arrival = hour(date_tz[is_arrival == T]) +
      minute(date_tz[is_arrival == T]) / 60 +
      second(date_tz[is_arrival == T]) / (60 * 60),
    .groups = "drop"
  ) %>%
  filter(!is.na(hour_day_departure) & !is.na(hour_day_arrival))

# compare the circular mean hour of departure between two types of trip
watson.williams.test(
  list(
    pb = hour_data %>%
      filter(trip == "Post Breeding") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours"),
    pm = hour_data %>%
      filter(trip == "Post Molting") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours")
  )
)

    Watson-Williams test for homogeneity of means

data:  pb and pm
F = 2.8203, df1 = 1, df2 = 214, p-value = 0.09454
sample estimates:
Circular Data: 
Type = angles 
Units = hours 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
mean of pb mean of pm 
 -2.819599  -4.515445 

Code

watson.wheeler.test(
  list(
    pb = hour_data %>%
      filter(trip == "Post Breeding") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours"),
    pm = hour_data %>%
      filter(trip == "Post Molting") %>%
      pull(hour_day_departure) %>%
      circular(., units = "hours")
  )
)

    Watson-Wheeler test for homogeneity of angles

data:  pb and pm
W = 1.1476, df = 2, p-value = 0.5634

More here: https://stats.stackexchange.com/questions/265865/comparisons-of-circular-means

4.2 Figure 1

Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.

Code

# to choose an individual ~180 dives for 36 hours
data_dive %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(DateTime), DateTime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(DateTime), DateTime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  )) %>%
  group_by(id, trip) %>%
  summarise(nb_dives = length(DiveNumber)) %>%
  arrange(nb_dives) %>%
  View()

Code

# import the csv export
tdr_data <- lapply(
  list.files(
    path = "../data/",
    pattern = "*_iknos_raw_data.csv",
    full.names = T
  ),
  function(x) {
    # extract id
    id_name <- last(strsplit(strsplit(x, split = "_")[[1]][1],
      split = "/"
    )[[1]])
    # load file
    dt <- readr::read_csv(x)
    # add id
    dt %>%
      mutate(id = as.numeric(id_name))
  }
) %>% bind_rows()

# choose the lag
nb_days <- 1.5

# compute the data
data_plot <- tdr_data %>%
  # convert time into date, a Posix
  mutate(
    DateTime = matlab_to_posix(time),
    date = as.Date(DateTime)
  ) %>%
  # add beginning and end date of the trip
  merge(
    .,
    data_dive %>%
      group_by(id, trip) %>%
      summarise(
        date_low = min(DateTime),
        date_high = max(DateTime),
        .groups = "drop"
      ),
    by = "id",
    all.x = TRUE
  ) %>%
  # identify what to keep
  mutate(to_keep = if_else(between(DateTime, date_low, date_high), T, F)) %>%
  # keep only between the beginning and the end
  filter(to_keep) %>%
  # sort to make sure the difftime are being done correctly
  arrange(id, time) %>%
  # merge to get lat and lon
  merge(
    .,
    data_dive %>%
      group_by(id, date) %>%
      summarise(
        lat = median(Lat, na.rm = T),
        lon = median(Lon, na.rm = T),
        .groups = "drop"
      ),
    by = c("id", "date"),
    all.x = TRUE
  ) %>%
  # then by individual
  group_by(id) %>%
  # calculate
  mutate( # the difference time with the first dive
    diff_start = abs(as.numeric(difftime(
      first(DateTime), DateTime,
      units = "days"
    ))),
    # the difference time with the last dive
    diff_end = abs(as.numeric(difftime(
      last(DateTime), DateTime,
      units = "days"
    )))
  ) %>%
  filter(diff_start <= nb_days | diff_end <= nb_days) %>%
  mutate(phase = factor(
    if_else(diff_start <= nb_days, "Departure", "Arrival"),
    levels = c("Departure", "Arrival")
  ))

# add bathy
setDT(data_plot)
setDT(data_dive)
setkeyv(data_plot, c("id", "DateTime"))
setkeyv(data_dive, c("id", "DateTime"))
data_plot <- copy(data_dive) %>%
  .[, .(id, DateTime, bathy, DiveTypeName)] %>%
  .[data_plot, roll = T]

# get sunrise sunset information
sun_data <- data_plot %>%
  select(date, lat, lon, id, phase) %>%
  unique() %>%
  ungroup() %>%
  mutate(getSunlightTimes(data = ., keep = c("sunrise", "sunset"))) %>%
  # fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1
  mutate(sunset = if_else( # if sunset is the next day
    abs(as.numeric(
      difftime(as.POSIXct(date), sunset, units = "days")
    )) > 1,
    # subtract 1
    sunset - (1 * 60 * 60 * 24),
    # otherwise leave it that way
    sunset
  ))

# trim if night beging before the trip or end after they arrive
sun_data <-
  merge(
    sun_data %>% select(date, lat, lon, id, phase, sunrise, sunset),
    data_plot %>%
      group_by(id, phase, trip) %>%
      summarise(
        date_low = min(DateTime),
        date_high = max(DateTime),
        .groups = "drop"
      ),
    by = c("id", "phase"),
    all.x = TRUE,
  ) %>%
  mutate(
    # make sure sunset is not outside the 36h period
    correct_sunset = if_else(
      # if not
      between(sunset, date_low, date_high),
      sunset,
      # if so, replace by one extremum
      if_else(sunset < date_low, date_low, date_high)
    ),
    # make sure sunrise is not outside the 36h period
    correct_sunrise = if_else(
      # if not
      between(sunrise, date_low, date_high),
      sunrise,
      # if so, replace by one extremum
      if_else(sunrise < date_low, date_low, date_high)
    )
  ) %>%
  # shift everything to start 2000-01-01 00:00:00
  mutate(
    correct_sunset_shift = correct_sunset - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    ),
    correct_sunrise_shift = correct_sunrise - (
      date_low - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
    )
  )

# plot
departure_dp <- data_plot %>%
  filter(phase == "Departure") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(DateTime = DateTime - (
    first(DateTime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = DateTime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Departure"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = DateTime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("+6h", "+12h", "+18h", "+24h", "+30h"),
    position = "top"
  ) +
  scale_y_reverse(limits = rev(range(data_plot$depth))) +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free_x",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    strip.text.y = element_blank(),
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      # arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_line(
      color = "black",
      arrow = arrow(
        length = unit(0.2, "lines"),
        type = "closed",
        ends = "first"
      )
    ),
    axis.title.x = element_blank()
  )

arrival_dp <- data_plot %>%
  filter(phase == "Arrival") %>%
  group_by(id, phase) %>%
  # make the date starts on 2000-01-01 00:00:00
  mutate(DateTime = DateTime - (
    first(DateTime) - as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
  )) %>%
  ungroup() %>%
  mutate(`Dive Type` = if_else(DiveTypeName == "Benthic", "Benthic", "Other")) %>%
  ggplot(aes(x = DateTime, y = depth, col = `Dive Type`, group = id)) +
  scale_color_manual(values = c("Benthic" = "#008c49", "Other" = "grey10")) +
  geom_rect(
    data = sun_data %>%
      filter(phase == "Arrival"),
    aes(
      xmin = correct_sunset_shift,
      xmax = correct_sunrise_shift,
      ymin = -Inf,
      ymax = +Inf
    ),
    alpha = 0.4,
    inherit.aes = FALSE
  ) +
  # geom_path(aes(x = DateTime, y = abs(bathy))) +
  geom_path() +
  scale_x_datetime(
    # date_breaks = "6 hours",
    # date_labels = "+ %Hh",
    breaks = seq.POSIXt(
      as.POSIXct("2000-01-01 06:00:00", tz = "UTC"),
      as.POSIXct("2000-01-02 06:00:00", tz = "UTC"),
      by = "6 hour"
    ),
    labels = c("-30h", "-24h", "-18h", "-12h", "-6h"),
    position = "top"
  ) +
  scale_y_reverse(limits = rev(range(data_plot$depth))) +
  labs(x = "Time", y = "Depth (m)") +
  # facet_nested(trip + id ~ phase,
  facet_nested(id ~ phase,
    scales = "free_x",
    # independent = "x"
  ) +
  theme(
    strip.placement = "outside",
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(color = "grey90"),
    axis.line.x = element_line(
      color = "black",
      arrow = arrow(length = unit(0.2, "lines"), type = "closed")
    ),
    axis.line.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

fig_label_x <-
  ggplot(data.frame(l = "Time", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

Code

guide_area() / fig_label_x / (departure_dp | arrival_dp) +
  plot_layout(
    guides = "collect",
    heights = c(0.5, 0.7, 6)
  ) &
  theme(
    legend.position = "top",
    plot.margin = unit(c(0, 0, 0, 0), "npc")
  )

Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.

Notes

• TODO: Add three dots between Departure and Arrival

4.3 Figure 2

Let’s create a table specific for this figure that only contains for each individual:

• the first dive
• the last dive
• all benthic dives

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # first we need to identify the start and the end of each trip by id
  group_by(id, trip) %>%
  mutate(
    is_arrival = if_else(DiveNumber == max(DiveNumber), T, F),
    is_departure = if_else(DiveNumber == min(DiveNumber), T, F)
  ) %>%
  # keep only the data of the first, last dives
  filter(is_arrival == T | is_departure == T) %>%
  # ungroup
  ungroup() %>%
  # then, since we're looking at local time, we keep data with location
  filter(!is.na(Lat)) %>%
  # change colnames to make getSunlightTimes works
  mutate(
    date = as.Date(date_tz),
    lat = Lat,
    lon = Lon
  ) %>%
  # calculate hour of sunset and hour of sunrise
  mutate(getSunlightTimes(
    data = .,
    keep = c("sunrise", "sunset"),
    tz = tz_lookup_coords(
      lat = Lat,
      lon = Lon,
      method = "accurate"
    )
  )) %>%
  # extract hours of sunset an sunrise
  mutate(
    h_sunrise = round(hour(sunrise) +
      minute(sunrise) / 60 +
      second(sunrise) / (60 * 60)),
    h_sunset = round(hour(sunset) +
      minute(sunset) / 60 +
      second(sunset) / (60 * 60))
  )

Code

# interesting circular time series
# https://campus.datacamp.com/courses/fraud-detection-in-r/introduction-motivation?ex=4

# for departure
windrose_departure_pb <- data_windrose %>%
  filter(is_departure == T & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n_distinct(id), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10),
    facet_wrap = "origin",
    hour_sunrise = data_windrose %>%
      filter(is_departure == T &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      psych::circadian.mean(.),
    hour_sunset = data_windrose %>%
      filter(is_departure == T &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      psych::circadian.mean(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  # geom_vline(xintercept = 6,
  #            linetype = "dashed") +
  geom_text(
    data = data.frame(
      x = rep(8, 2),
      y = c(6, 11),
      text = c(5, 10)
    ),
    aes(x = x, y = y, label = text),
    size = 2
  )

windrose_departure_pm <- data_windrose %>%
  filter(is_departure == T & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n_distinct(id), .groups = "drop") %>%
  mutate(origin = "Departure") %>%
  windrose(.,
    grid = c(5, 10),
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(is_departure == T &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      psych::circadian.mean(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(is_departure == T &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      psych::circadian.mean(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
# geom_vline(xintercept = 6,
#            linetype = "dashed") +
# geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

# for arrival
windrose_arrival_pb <- data_windrose %>%
  filter(is_arrival == T & trip == "Post Breeding") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n_distinct(id), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10),
    legend_position = "top",
    facet_grid = T,
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(is_arrival == T &
        trip == "Post Breeding") %>%
      pull(h_sunrise) %>%
      psych::circadian.mean(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(is_arrival == T &
        trip == "Post Breeding") %>%
      pull(h_sunset) %>%
      psych::circadian.mean(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

windrose_arrival_pm <- data_windrose %>%
  filter(is_arrival == T & trip == "Post Molting") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time, trip) %>%
  summarise(nb_ind_hour = n_distinct(id), .groups = "drop") %>%
  mutate(origin = "Arrival") %>%
  windrose(.,
    grid = c(5, 10),
    legend_position = "top",
    facet_wrap = "trip",
    hour_sunrise = data_windrose %>%
      group_by(id) %>%
      filter(is_arrival == T &
        trip == "Post Molting") %>%
      pull(h_sunrise) %>%
      psych::circadian.mean(.),
    hour_sunset = data_windrose %>%
      group_by(id) %>%
      filter(is_arrival == T &
        trip == "Post Molting") %>%
      pull(h_sunset) %>%
      psych::circadian.mean(.)
  ) +
  theme(
    axis.title = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) #+
#   geom_text(data = data.frame(
#   x = rep(8,3),
#   y = c(6, 11, 16),
#   text = c(5, 10, 15)
# ),
# aes(x = x, y = y, label = text),
# size = 2)

label_x <- ggplot(data.frame(l = "Time (hours)", x = 1, y = 1)) +
  geom_text(aes(x, y, label = l)) +
  theme_void() +
  coord_cartesian(clip = "off")

# combine all
(windrose_departure_pb + windrose_arrival_pb) /
  (windrose_departure_pm + windrose_arrival_pm) /
  label_x +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  plot_layout(
    heights = c(2, 2, 0),
    guides = "collect"
  ) &
  # plot_annotation(tag_levels = c("A")) &
  # theme(plot.tag.position = c(-0.03, 0.8),
  #       plot.tag = element_text(size = 8, hjust = 0, vjust = 0)) +
  theme(legend.position = "top")

Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.

Note

• n = 361, this is the number of individuals having location data (allowing to calculate the local time) for their first and/or their last dives

Code

# departure PM
pm_dep_mean <- psych::circadian.mean(
  data_windrose %>%
    filter(is_departure == T &
      trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# departure PB
pb_dep_mean <- psych::circadian.mean(
  data_windrose %>%
    filter(is_departure == T &
      trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# arrival PM
pm_arr_mean <- psych::circadian.mean(
  data_windrose %>%
    filter(is_arrival == T &
      trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

# arrival PB
pb_arr_mean <- psych::circadian.mean(
  data_windrose %>%
    filter(is_arrival == T &
      trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>% pull(time)
)

Code

# functio to convert time in radian (probably already exist in some package...)
hour_time_2_radian <- function(x, high = 24, low = 0) {
  # https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles
  ptp <- high - low
  rad <- x * (2 * pi) / ptp
  rad <- (rad + pi) %% (2 * pi) - pi
  return(rad)
}

Code

# departure
test_departure_pm <- CircStats::r.test(
  data_windrose %>%
    filter(is_departure == T &
      trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)
test_departure_pb <- CircStats::r.test(
  data_windrose %>%
    filter(is_departure == T &
      trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)

# arrival
test_arrival_pm <- CircStats::r.test(
  data_windrose %>%
    filter(is_arrival == T &
      trip == "Post Molting") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)
test_arrival_pb <- CircStats::r.test(
  data_windrose %>%
    filter(is_arrival == T &
      trip == "Post Breeding") %>%
    mutate(time = hour(date_tz) +
      minute(date_tz) / 60 +
      second(date_tz) / (60 * 60)) %>%
    pull(time) %>%
    hour_time_2_radian(.)
)

Code

data.frame(
  Departure = c(
    paste0(
      round(pb_dep_mean, 1),
      " (R = ",
      round(test_departure_pb$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_departure_pb$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_departure_pb$p.value, 2))
      ),
      ")"
    ),
    paste0(
      round(pm_dep_mean, 1),
      " (R = ",
      round(test_departure_pm$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_departure_pm$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_departure_pm$p.value, 2))
      ),
      ")"
    )
  ),
  Arrival = c(
    paste0(
      round(pb_arr_mean, 1),
      " (R = ",
      round(test_arrival_pb$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_arrival_pb$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_arrival_pb$p.value, 2))
      ),
      ")"
    ),
    paste0(
      round(pm_arr_mean, 1),
      " (R = ",
      round(test_arrival_pm$r.bar, 2),
      ", p-value ",
      ifelse(
        round(test_arrival_pm$p.value, 2) == 0,
        "< 0.001",
        paste("= ", round(test_arrival_pm$p.value, 2))
      ),
      ")"
    )
  ),
  rownames = c("Post Breeding", "Post Molting")
) %>%
  gt(rowname_col = "rownames") %>%
  cols_label(
    Arrival = "Hour of Arrival",
    Departure = "Hour of Departure"
  )

Table 3:

(Not meant to be in the paper) Average hour of departure and arrival for post breeding and post molting trip with the result of the associated Rayleigh test.

	Hour of Departure	Hour of Arrival
Post Breeding	20.6 (R = 0.22, p-value < 0.001)	9.3 (R = 0.02, p-value = 0.93)
Post Molting	19.3 (R = 0.23, p-value < 0.001)	13.7 (R = 0.15, p-value = 0.04)

Good ressource to interpret Rayleigh test: https://stats.stackexchange.com/questions/433396/interpretation-of-rayleigh-test-results

4.4 Figure 3

Code

# with or without bathymetry
with_bathy <- T

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_in.rds")) {
    trip_zoom_in <- readRDS("../export/background_ggoceanmap_zoom_in.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_in <- basemap(
      # limits = c(170, -110, 30, 59),
      # limits = c(-155, -110, 30, 60),
      limits = c(-165, -110, 30, 65),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_in)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_in)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_in <- trip_zoom_in +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_in$layers <- trip_zoom_in$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_in, "../export/background_ggoceanmap_zoom_in.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_in <- basemap(
    # limits = c(170, -110, 20, 59),
    # limits = c(-155, -110, 30, 60),
    limits = c(-165, -110, 30, 65),
    land.col = "black",
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_in)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_in)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_in <- trip_zoom_in +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_in$layers <- trip_zoom_in$layers[c(2, 1)]
}

Code

# data_dive_longer <- pivot_longer(
#   data_dive %>%
#     # keep data with locations
#     filter(!is.na(Lat)),
#   c(
#     "shark_area",
#     "orca_area",
#     "overlap_area",
#     "total_area"
#   ),
#   names_to = "areas",
#   values_to = "in_out"
# )
#
# data_descriptive_ind <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, id) %>%
#   summarise(
#     n_tot = n(),
#     n_inside = length(id[in_out == 1]),
#     n_outside = length(id[in_out == 0]),
#     .groups = "drop"
#   ) %>%
#   mutate(
#     perc_inside = n_inside / n_tot,
#     perc_outside = n_outside / n_tot
#   )
#
# data_descriptive_area <- data_descriptive_ind %>%
#   group_by(areas) %>%
#   summarise(
#     mean_perc_inside = wtd.mean(perc_inside, n_tot),
#     mean_perc_outside = wtd.mean(perc_outside, n_tot),
#     sd_perc_inside = sqrt(wtd.var(perc_inside, n_tot)),
#     sd_perc_outside = sqrt(wtd.var(perc_outside, n_tot)),
#     # based from here: https://stats.stackexchange.com/questions/25895/computing-standard-error-in-weighted-mean-estimation
#     se_perc_inside = sqrt(wtd.var(perc_inside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#     se_perc_outside = sqrt(wtd.var(perc_outside, n_tot)*sum((n_tot/sum(n_tot))^2)),
#   ) %>%
#   pivot_longer(
#     cols = ends_with("inside") |
#       ends_with("outside"),
#     names_to = c("stat", "in_out"),
#     names_pattern = "(.*)_perc_(.*)"
#   ) %>%
#   pivot_wider(
#     id_cols = c("areas", "in_out"),
#     names_from = "stat",
#     values_from = "value"
#   ) %>%
#   mutate(areas = factor(areas,
#                         level = c(
#                           "overlap_area",
#                           "orca_area",
#                           "shark_area",
#                           "total_area"
#                         )
#   ))
#
# data_stat_area <- data_dive_longer %>%
#   filter(DiveTypeName == "Benthic") %>%
#   group_by(areas, in_out) %>%
#   summarise(nb_area = n(), .groups = "drop_last") %>%
#   mutate(nb_total = sum(nb_area)) %>%
#   summarise(rstatix::prop_test(x = nb_area, n = nb_total))

Code

# inside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside", -mean, mean),
#     sd = if_else(in_out == "outside", -sd, sd),
#     se = if_else(in_out == "outside", -se, se)
#   ) %>%
#   filter(in_out == "inside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(
#     show.legend = F,
#     mapping = aes(
#       x = mean,
#       y = areas,
#       fill = in_out),
#     alpha = 0.3
#   ) +
#   scale_fill_manual(values = c("grey70")) +
#   geom_errorbar(mapping = aes(xmin = -se,
#                               xmax = se),
#                 # reduce size of horizontal bar of the errorbar
#                 width = 0.05) +
#   coord_flip(xlim = c(0, max(data_descriptive_area$mean))) +
#   scale_x_continuous(
#     # breaks = c(-0.6,-0.4,-0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0,0.1), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x)
#         (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   theme(
#     panel.spacing.y = unit(0, "mm"),
#     panel.grid.major = element_blank(),
#     panel.grid.minor = element_blank(),
#     panel.background = element_blank(),
#     axis.line.y = element_line(arrow = arrow(
#       angle = 20,
#       length = unit(.10, "inches"),
#       type = "closed"
#     )),
#     axis.title = element_blank(),
#     axis.ticks.length.x = unit(0, units = "mm"),
#     axis.text.x = element_blank()
#   )
#
# outside_bar_plot = data_descriptive_area %>%
#   mutate(
#     mean = if_else(in_out == "outside",-mean, mean),
#     sd = if_else(in_out == "outside",-sd, sd),
#     se = if_else(in_out == "outside",-se, se)
#   ) %>%
#   filter(in_out == "outside") %>%
#   ggplot(aes(x = mean, y = areas, fill = areas)) +
#   facet_grid(in_out ~ ., scales = "free_y") +
#   geom_col(show.legend = F) +
#   scale_fill_manual(values = c("#9A6BEC", "#E66100", "#E1E826", "#004D40")) +
#   new_scale_fill() +
#   geom_col(show.legend = F, mapping = aes(x = mean,
#                                           y = areas,
#                                           fill = in_out),
#                                           alpha = 0.3) +
#   scale_fill_manual(values = c("grey30")) +
#   geom_errorbar(
#     mapping = aes(
#       xmin = -se,
#       xmax = se
#     ),
#     # reduce size of horizontal bar of the errorbar
#     width = 0.05
#   ) +
#   coord_flip(xlim = c(-max(data_descriptive_area$mean), 0)) +
#   scale_x_continuous(
#     # breaks = c(-0.6, -0.4, -0.2, 0.2, 0.4, 0.6),
#     expand = expansion(mult = c(0.1,0), add = 0),
#     labels = function(x) {
#       scales::percent(abs(x), 1)
#     }
#   ) +
#   scale_y_discrete(
#     labels = function(x) {
#       lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1])
#     }
#   ) +
#   labs(y = "Areas") +
#   theme(panel.spacing.y = unit(0, "mm"),
#         panel.grid.major = element_blank(),
#         panel.grid.minor = element_blank(),
#         panel.background = element_blank(),
#         axis.line.y = element_line(arrow = arrow(
#           angle = 20,
#           length = unit(.10, "inches"),
#           type = "closed",
#           ends = "first"
#         )),
#         axis.ticks.length.x = unit(0,units = "mm"),
#         axis.title.y = element_blank())
#
# label_x_bar_plot = ggplot(data.frame(l = "Proportion of benthic dives",
#                                      x = 1,
#                                      y = 1)) +
#   geom_text(aes(x, y, label = l, angle = 90)) +
#   theme(aspect.ratio = 10,
#         axis.title = element_blank(),
#         axis.text = element_blank(),
#         axis.ticks = element_blank()) +
#   theme_void()

Code

# color palette
# https://davidmathlogic.com/colorblind/#%23D81B60-%231E88E5-%23FFC107-%23004D40
# combine areas
areas_coord <- rbind(
  Orca_area_coord_to_display %>% mutate(Predator = "orca"),
  Shark_area_coord_to_display %>% mutate(Predator = "shark"),
  Overlap_area_coord_to_display %>% mutate(Predator = "overlap")
) %>%
  mutate(Predator = factor(Predator, levels = c("orca", "shark", "overlap")))

# data use to compute kernel density estimation
df_kernel_dens_benthic <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Lon, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

# transform data into sf object
df_kernel_dens_sf_benthic <- st_as_sf(
  df_kernel_dens_benthic,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf_benthic <- group_by(df_kernel_dens_sf_benthic, origin)

# kernel density estimation
df_kernel_dens_sf_kde_benthic <- st_kde(
  df_kernel_dens_sf_benthic,
  # H = diag(c(
  #   MASS::bandwidth.nrd(
  #     sf::st_coordinates(df_kernel_dens_sf_benthic)[, 1]
  #   ),
  #   MASS::bandwidth.nrd(
  #     sf::st_coordinates(df_kernel_dens_sf_benthic)[, 2]
  #   )
  # ) / 4)^2
)

# df_kernel_dens_sf_kde_benthic <- st_kde(df_kernel_dens_sf_benthic,
#                                         H = ks::Hpi(sf::st_coordinates(df_kernel_dens_sf_benthic),
#                                                          pilot = "dunconstr"))

# df_kernel_dens_sf_kde_benthic <- st_kde(df_kernel_dens_sf_benthic,
#                                         H = ks::Hlscv.diag(sf::st_coordinates(df_kernel_dens_sf_benthic)))

# df_kernel_dens_sf_kde_benthic <- st_kde(df_kernel_dens_sf_benthic,
#                                 H = diag(c(
#                                   MASS::bandwidth.nrd(
#                                     sf::st_coordinates(df_kernel_dens_sf_benthic)[, 1]
#                                   ),
#                                   MASS::bandwidth.nrd(
#                                     sf::st_coordinates(df_kernel_dens_sf_benthic)[, 2]
#                                   )
#                                 ))
# )


# https://loading.io/color/feature/BuGn-9/
cols <- c(
  "99%" = "#ccece6",
  "95%" = "#99d8c9",
  "80%" = "#41ae76",
  "50%" = "#00441b"
)

p1 <- trip_zoom_in +
  # remove legend
  guides(fill = "none") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # legend
  labs(fill = "Probability contours") +
  # legend at the top
  theme(legend.position = "top") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# Extract the size legend - leg1
leg1 <- gtable_filter(ggplot_gtable(ggplot_build(p1)), "guide-box")

p2 <- trip_zoom_in +
  # tweak legend
  guides(fill = guide_legend(
    override.aes = list(
      colour = "grey60"
    )
  )) +
  theme(
    legend.key = element_rect(color = "grey60"),
    legend.text = element_text(colour = "grey20"),
    legend.title = element_text(colour = "grey20"),
    legend.justification = c("right", "top"),
    legend.position = c(0.99, 0.99),
    legend.background = element_rect(fill = alpha("white", 0))
  ) +
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour( # geospatial kernel
      df_kernel_dens_sf_kde_benthic,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  # scale_fill_viridis_d(option = "magma", direction = -1) +
  # scale_fill_brewer(palette = "heat")
  # scale_fill_brewer(palette = "BuGn") +
  scale_fill_manual(values = cols) +
  # scale_fill_grey() +
  # remove legend
  guides(fill = "none") +
  # # no display alpha
  # guides(alpha = "none") +
  # new scale
  new_scale_fill() +
  # add spatial polygon
  geom_spatial_polygon(
    data = areas_coord,
    aes(
      x = lon,
      y = lat,
      fill = Predator
    ),
    alpha = 0.35,
    crs = 4326
  ) +
  # color palette
  scale_fill_manual(values = c("darkred", "#E1E826", "#E66100")) +
  scale_color_manual(values = c("darkred", "#E1E826", "#E66100")) +
  # remove legend for areas
  guides(fill = "none", col = "none") +
  # set limit from predator_map$layers[[1]]$data
  coord_sf(
    xlim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$xmax
    ),
    ylim = c(
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymin,
      st_bbox(trip_zoom_in$layers[[1]]$data)$ymax
    )
  )


# reordering layers (continent on top)
p2$layers <- if (with_bathy) {
  p2$layers[c(1, 2, 5, 4, 3)]
} else {
  p2$layers[c(1, 4, 3, 2)]
}

# convert into ggplot
predator_map <- as_ggplot(leg1) / p2 +
  plot_layout(
    ncol = 1,
    heights = c(1, 12)
  )

Code

# load image icon
shark_img <- readPNG("../img/shark_custom.png", native = T)
orca_img <- readPNG("../img/orca_custom.png", native = T)

# draw the map
predator_map +
  inset_element(
    p = orca_img,
    left = 0.2,
    bottom = 0.82,
    right = 0.3,
    top = 0.92
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  ) +
  inset_element(
    p = shark_img,
    left = 0.15,
    bottom = 0.15,
    right = 0.25,
    top = 0.25
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  )

Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.

Code

# draw the map
predator_map +
  ggspatial::geom_spatial_point(
    data = data_dive %>%
      filter(!is.na(Lat) &
        DiveTypeName == "Benthic"),
    aes(x = Lon, y = Lat), alpha = 0.2, size = 1
  ) +
  guides(alpha = "none", size = "none") +
  inset_element(
    p = orca_img,
    left = 0.2,
    bottom = 0.82,
    right = 0.3,
    top = 0.92
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  ) +
  inset_element(
    p = shark_img,
    left = 0.15,
    bottom = 0.15,
    right = 0.25,
    top = 0.25
  ) +
  theme(
    rect = element_rect(fill = "transparent"),
    plot.background = element_rect(color = NA)
  )

Figure 4: Same as the previous one, but with the benthic dives (black dots). We can see that some benthic happen outside orca areas…

Notes

• ?
• n = 392, this is the number of individuals having at least one benthic with one location associated with (filter(!is.na(Lat) & DiveTypeName == "Benthic"))

4.5 Figure 4

Code

# https://mycolor.space/?hex=%23008C49&sub=1
# build palette color
cols_ind_1 <- c("#008c49", "#00c9d3", "#bffcfd", "#458081")
# test for color blind
# => https://davidmathlogic.com/colorblind/#%23008C49-%2300C9D3-%23BFFCFD-%23458081
cols_ind_2 <- c("#008c49", "#002e00", "#00c9d3", "#deb887")
# test for color blind
# => https://davidmathlogic.com/colorblind/#%23008C49-%23002E00-%2300C9D3-%23DEB887

Code

# main plot
bathy_prop <- data_dive %>%
  # keep only rows with location data
  filter(!is.na(Lat)) %>%
  # keep only negative bathy
  filter(bathy < 0) %>%
  # and remove outliers
  filter(bathy > -6000) %>%
  # positive bathy
  mutate(bathy = abs(bathy)) %>%
  # create class of bathymetry
  mutate(bathy_class = cut(
    bathy,
    seq(0, 6000, 400),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(bathy_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per bathy_class
  group_by(bathy_class) %>%
  # to get the Proportion of different dive types per bathy_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = bathy_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Bathymetry  (m)",
    y = "Dive type proportion (%)"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "top",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # second plot
# bathy_hist <- data_dive %>%
#   # keep only negative bathy
#   filter(bathy < 0) %>%
#   # and remove outliers
#   filter(bathy > -6000) %>%
#   # create class of bathymetry
#   mutate(bathy_class = fct_rev(cut(
#     bathy,
#     seq(-6000, 0, 400),
#     ordered_result = T,
#     dig.lab = 4
#   ))) %>%
#   # calculate by bath_class and animal
#   group_by(bathy_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = bathy_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))

Code

# main plot
dist_coast_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_coast > 0) %>%
  # remove outliers
  filter(dist_coast * 1000 * 111 <= 1900) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    # convert decimal degree/1000
    dist_coast * 1000 * 111,
    seq(0, 1900, 100),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from the coast (km)",
    y = "Dive type proportion"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # the second plot
# dist_coast_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_coast * 1000 * 111 <= 1900) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     # convert decimal degree/1000
#     dist_coast * 1000 * 111,
#     seq(0, 2000, 100),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# the actual plot
# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))

Code

# the main plot
dist_dep_prop <- data_dive %>%
  # filter out animal without dist_coast
  filter(dist_dep > 0) %>%
  # remove outliers
  filter(dist_dep / 1000 < 4200) %>%
  # create class of bathymetry
  mutate(dist_class = cut(
    dist_dep / 1000,
    seq(0, 4200, 300),
    ordered_result = T,
    dig.lab = 4
  )) %>%
  # calculate by bath_class and animal
  group_by(dist_class, DiveTypeName) %>%
  # the number of dives
  summarise(N = n(), .groups = "drop") %>%
  # and divide by the total number of dives per dist_class
  group_by(dist_class) %>%
  # to get the Proportion of different dive types per dist_class
  mutate(Proportion = N / sum(N)) %>%
  # ungroup => not required but that let the dataset clean
  ungroup() %>%
  # then plot
  ggplot(aes(x = dist_class, y = Proportion)) +
  # the area
  geom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +
  # orientation of x labels
  guides(x = guide_axis(angle = 45)) +
  # format y axis
  scale_y_continuous(
    labels = function(x) {
      paste0(x * 100, "%")
    }
  ) +
  labs(
    x = "Distance from A<U+00F1>o Nuevo (km)",
    y = "Dive type proportion"
  ) +
  # scale_fill_viridis(
  #   "Dive Type",
  #   option = "plasma",
  #   discrete = T,
  #   direction = -1
  # ) +
  scale_fill_manual("Dive Type",
    values = cols_ind_1
  ) +
  # position legend at the top
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 15)),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    )),
    plot.tag.position = c(-0.03, 1)
  )

# # the second plot
# dist_dep_hist <- data_dive %>%
#   # filter out animal without dist_coast
#   filter(dist_coast > 0) %>%
#   # remove outliers
#   filter(dist_dep / 1000 < 4200) %>%
#   # create class of bathymetry
#   mutate(dist_class = cut(
#     dist_dep / 1000,
#     seq(0, 4200, 300),
#     ordered_result = T,
#     dig.lab = 4
#   )) %>%
#   # calculate by bath_class and animal
#   group_by(dist_class, DiveTypeName) %>%
#   # the number of dives
#   summarise(N = n(), .groups = "drop") %>%
#   # ggplot
#   ggplot(aes(
#     x = dist_class,
#     y = N,
#     fill = DiveTypeName,
#     width = 0.5
#   )) +
#   geom_bar(
#     stat = "identity",
#     position = "dodge"
#   ) +
#   scale_fill_viridis(
#     option = "plasma",
#     discrete = T,
#     direction = -1
#   ) +
#   theme_void() +
#   theme(legend.position = "top")
#
# # the actual plot
# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))

Code

# y axis label position
ylab <- dist_coast_prop$labels$y
bathy_prop$labels$y <- dist_dep_prop$labels$y <- " "

Code

# display
(guide_area() + bathy_prop + dist_coast_prop + dist_dep_prop) +
  theme(plot.margin = unit(c(0, .02, -1, .05), units = "npc")) +
  plot_annotation(tag_levels = "A") +
  plot_layout(nrow = 5, guides = "collect")
# grid::grid.draw(grid::textGrob(ylab, x = 0.07, y = 0.55, rot = 90))

Figure 5: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).

Notes

• n = 406, because we removed individuals without location data (filter(!is.na(Lat)))

4.6 Figure 5

Code

# the last plot
prop_at_sea_dt <- data_dive %>%
  group_by(id) %>%
  # time difference between the first date and the others (telling R to take diff b/w first date and all other dates)
  mutate(nb_days_departure = trunc(as.numeric(difftime(
    date, first(date),
    units = "days"
  )))) %>%
  # group by day/dive_type/seal
  group_by(nb_days_departure, DiveTypeName, id) %>%
  # count number of dives
  summarise(nb_daily_dives_type = n(), .groups = "drop") %>%
  # group by day/seal
  group_by(nb_days_departure, id) %>%
  # count the total number of dives per day/seal
  mutate(nb_daily_dives = sum(nb_daily_dives_type)) %>%
  # order by seals/date
  arrange(id, nb_days_departure) %>%
  # perform our calculation per seal
  group_by(id) %>%
  # calculate of the Proportion of time since departure
  mutate(perc_time_at_sea = round(
    nb_days_departure / max(nb_days_departure) * 100, 1
  )) %>%
  # filter on benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # add class of Proportion of time at sea
  mutate(day_class = cut(perc_time_at_sea, seq(0, 100, 5),
    include.lowest = T
  )) %>%
  # by class of Proportion of time at sea
  group_by(day_class) %>%
  # calculate the proportion of daily benthic dives
  summarise(
    nb_benthic_class = sum(nb_daily_dives_type),
    nb_total_class = sum(nb_daily_dives),
    prop_class = nb_benthic_class / nb_total_class
  )

# plot
prop_at_sea <- prop_at_sea_dt %>%
  # intiate plot
  ggplot(
    .,
    aes(x = day_class, y = prop_class)
  ) +
  geom_bar(stat = "identity", fill = "#008c49") +
  guides(x = guide_axis(angle = 45)) +
  scale_y_continuous(
    limits = c(0, 1),
    # breaks=seq(0,0.6,0.1),
    labels = function(x) {
      scales::percent(abs(x), 1)
    }
  ) +
  annotate("text", x = 5.3, y = 0.4, label = "Departure") +
  annotate("text", x = 15, y = 0.43, label = "Arrival") +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(4, 17),
  #     x2 = c(1.5, 19.5),
  #     y1 = c(0.4, 0.5),
  #     y2 = c(0.31, 0.55)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = 0.3
  # ) +
  geom_curve(
    data = tibble(
      x1 = c(4),
      x2 = c(1),
      y1 = c(0.4),
      y2 = c(0.33)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = 0.3
  ) +
  geom_curve(
    data = tibble(
      x1 = c(16),
      x2 = c(19.7),
      y1 = c(0.43),
      y2 = c(0.5)
    ),
    aes(
      x = x1,
      y = y1,
      xend = x2,
      yend = y2
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    linewidth = 0.5,
    color = "gray20",
    curvature = -0.56
  ) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  ) +
  labs(
    x = "Time spent at sea (%)",
    y = "Proportion of benthic dives"
  )

Code

# display
prop_at_sea

Figure 6: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.

Notes

• n = 477 all individuals (since they all have at leat one benthic dives)

Test to compare if the proportion of benthic dives at arrival is higher than the one at departure:

• \(H_{0}: p_{departure} \ge p_{arrival}\)
• \(H_{1}: p_{departure} \lt p_{arrival}\)

Code

# test to compare the proportion between arrival and departure
z_test <- prop.test(
  x = c(
    prop_at_sea_dt %>%
      filter(day_class == first(day_class)) %>%
      pull(nb_benthic_class),
    prop_at_sea_dt %>%
      filter(day_class == last(day_class)) %>%
      pull(nb_benthic_class)
  ),
  n = c(
    prop_at_sea_dt %>%
      filter(day_class == first(day_class)) %>%
      pull(nb_total_class),
    prop_at_sea_dt %>%
      filter(day_class == last(day_class)) %>%
      pull(nb_total_class)
  ),
  alternative = "less"
)
z_test

    2-sample test for equality of proportions with continuity correction

data:  c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_benthic_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_benthic_class)) out of c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_total_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_total_class))
X-squared = 4305.3, df = 1, p-value < 2.2e-16
alternative hypothesis: less
95 percent confidence interval:
 -1.0000000 -0.1808821
sample estimates:
   prop 1    prop 2 
0.2880853 0.4734710 

Code

# extract the z statistic
sqrt(z_test$statistic)

X-squared 
 65.61442 

\(H_{0}\) is rejected, so we can say that the porportion of benthic dives at arrival is significantly superior than departure!

Based on: http://www.sthda.com/english/wiki/comparing-proportions-in-r

and https://bookdown.org/danieljcarter/r4steph/z-test-for-proportions.html

4.7 Supplementary

4.7.1 Table 1

Code

# shark and orca area = NA if no location data
data_dive <- data_dive %>%
  mutate(
    shark_area = if_else(is.na(Lat), NA, shark_area),
    orca_area = if_else(is.na(Lat), NA, orca_area)
  ) %>%
  # add year of deployment
  group_by(id) %>%
  mutate(
    year = first(format(date, format = "%Y")),
    shark_and_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area == 2, 1, 0)
    ),
    shark_or_orca = if_else(is.na(Lat),
      NA,
      if_else(shark_area + orca_area >= 1, 1, 0)
    )
  )

# build table
data_dive %>%
  group_by(year, id) %>%
  summarise(
    n_dives_ind_all = length(DiveNumber),
    n_dives_ind_w_loc = length(DiveNumber[!is.na(Lat)]),
    n_dives_ind_wo_loc = length(DiveNumber[is.na(Lat)]),
    n_benthic_dives_all = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic"])
    ),
    n_benthic_dives_w_loc = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" & !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_shark = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_in_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        orca_area == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_and_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_and_orca == 1 &
        !is.na(Lat)])
    ),
    n_benthic_dives_ind_shark_or_orca = if_else(length(id[DiveTypeName == "Benthic"]) == 0,
      NA_real_,
      length(DiveNumber[DiveTypeName == "Benthic" &
        shark_or_orca == 1 &
        !is.na(Lat)])
    ),
    hour_day_departure = hour(date_tz[DiveNumber == min(DiveNumber)]) +
      minute(date_tz[DiveNumber == min(DiveNumber)]) / 60 +
      second(date_tz[DiveNumber == min(DiveNumber)]) / (60 * 60),
    hour_day_arrival = hour(date_tz[DiveNumber == max(DiveNumber)]) +
      minute(date_tz[DiveNumber == max(DiveNumber)]) / 60 +
      second(date_tz[DiveNumber == max(DiveNumber)]) / (60 * 60),
    .groups = "drop"
  ) %>%
  group_by(year) %>%
  summarise(
    n_ind = length(id),
    n_dives_all = sum(n_dives_ind_all, na.rm = T),
    n_dives_w_loc = sum(n_dives_ind_w_loc, na.rm = T),
    n_dives_wo_loc = sum(n_dives_ind_wo_loc, na.rm = T),
    n_benthic_dives = sum(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_proportion = sum(n_benthic_dives, na.rm = T) / sum(n_dives_all, na.rm = T),
    n_benthic_dives_in_shark_proportion = sum(n_benthic_dives_ind_in_shark, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_orca_proportion = sum(n_benthic_dives_ind_in_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_and_orca = sum(n_benthic_dives_ind_shark_and_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_benthic_dives_in_shark_or_orca = sum(n_benthic_dives_ind_shark_or_orca, na.rm = T) / sum(n_benthic_dives_w_loc, na.rm = T),
    n_dives_mean = mean(n_dives_ind_all, na.rm = T),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = sd(n_dives_ind_all, na.rm = T),
    n_benthic_dives_mean = mean(n_benthic_dives_all, na.rm = T),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = sd(n_benthic_dives_all, na.rm = T),
    hour_day_departure_mean = psych::circadian.mean(hour_day_departure, na.rm = T),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = psych::circadian.sd(hour_day_departure, na.rm = T)$sd,
    hour_day_arrival_mean = psych::circadian.mean(hour_day_arrival, na.rm = T),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = psych::circadian.sd(hour_day_arrival, na.rm = T)$sd
  ) %>%
  gt() %>%
  tab_spanner(
    label = "All dives",
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc"
    )
  ) %>%
  tab_spanner(
    label = "Benthic dives",
    columns = c(
      "n_benthic_dives",
      "n_benthic_dives_proportion"
    )
  ) %>%
  tab_spanner(
    label = "% Benthic dives",
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    )
  ) %>%
  tab_spanner(
    label = "# dives",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "# benthic dives",
    columns = c(
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour departure",
    columns = c(
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd"
    )
  ) %>%
  tab_spanner(
    label = "Hour arrival",
    columns = c(
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  tab_spanner(
    label = "per individual",
    columns = c(
      "n_dives_mean",
      "n_dives_plus_minus",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_plus_minus",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_plus_minus",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_plus_minus",
      "hour_day_arrival_sd"
    )
  ) %>%
  # rename columns
  cols_label(
    year = "Year",
    n_ind = "# seals",
    n_dives_all = "Total",
    n_dives_w_loc = "w/ loc",
    n_dives_wo_loc = "w/o loc",
    n_benthic_dives = "Total",
    n_benthic_dives_proportion = "Proportion",
    n_benthic_dives_in_shark_proportion = "Shark area",
    n_benthic_dives_in_orca_proportion = "Orca area",
    n_benthic_dives_in_shark_or_orca = "Total",
    n_benthic_dives_in_shark_and_orca = "Overlap",
    n_dives_mean = md("Mean"),
    n_dives_plus_minus = "\u00b1",
    n_dives_sd = md("SD"),
    n_benthic_dives_mean = md("Mean"),
    n_benthic_dives_plus_minus = "\u00b1",
    n_benthic_dives_sd = md("SD"),
    hour_day_departure_mean = md("Mean"),
    hour_day_departure_plus_minus = "\u00b1",
    hour_day_departure_sd = md("SD"),
    hour_day_arrival_mean = md("Mean"),
    hour_day_arrival_plus_minus = "\u00b1",
    hour_day_arrival_sd = md("SD")
  ) %>%
  fmt_percent(
    columns = c(
      "n_benthic_dives_in_shark_proportion",
      "n_benthic_dives_in_orca_proportion",
      "n_benthic_dives_proportion",
      "n_benthic_dives_in_shark_or_orca",
      "n_benthic_dives_in_shark_and_orca"
    ),
    decimal = 1
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_all",
      "n_dives_w_loc",
      "n_dives_wo_loc",
      "n_benthic_dives"
    ),
    decimal = 0
  ) %>%
  fmt_number(
    columns = c(
      "n_dives_mean",
      "n_dives_sd",
      "n_benthic_dives_mean",
      "n_benthic_dives_sd",
      "hour_day_departure_mean",
      "hour_day_departure_sd",
      "hour_day_arrival_mean",
      "hour_day_arrival_sd"
    ),
    decimal = 1,
    use_seps = FALSE
  ) %>%
  # alignement
  cols_align(
    columns = c(
      "n_dives_sd",
      "n_benthic_dives_sd",
      "hour_day_departure_sd",
      "hour_day_arrival_sd"
    ),
    align = "left"
  ) %>%
  # color rows
  opt_row_striping() %>%
  tab_options(table.width = pct(180))

											per individual
Year	# seals	All dives			Benthic dives		% Benthic dives				# dives			# benthic dives			Hour departure			Hour arrival
		Total	w/ loc	w/o loc	Total	Proportion	Shark area	Orca area	Total	Overlap	Mean	±	SD	Mean	±	SD	Mean	±	SD	Mean	±	SD
2004	23	244,042	157,085	86,957	7,754	3.2%	47.7%	93.1%	93.3%	47.5%	10610.5	±	3052.7	337.1	±	393.3	16.6	±	1.6	16.6	±	1.2
2005	38	311,574	267,603	43,971	14,438	4.6%	57.1%	66.3%	66.6%	56.7%	8199.3	±	3621.0	379.9	±	847.9	19.5	±	1.5	19.1	±	1.8
2006	34	277,497	227,242	50,255	5,296	1.9%	63.0%	94.5%	94.5%	62.9%	8161.7	±	3706.1	155.8	±	205.0	23.7	±	1.8	12.4	±	2.4
2007	36	271,346	165,371	105,975	7,643	2.8%	79.2%	81.3%	82.3%	78.2%	7537.4	±	4077.9	212.3	±	401.1	17.9	±	1.6	0.0	±	2.3
2008	35	250,958	181,586	69,372	4,636	1.8%	82.3%	86.3%	86.4%	82.2%	7170.2	±	3508.0	132.5	±	176.6	12.6	±	2.0	12.0	±	1.8
2009	25	150,609	92,402	58,207	4,280	2.8%	73.9%	85.6%	85.6%	73.9%	6024.4	±	3172.6	171.2	±	197.4	18.7	±	1.4	16.0	±	1.7
2010	38	263,334	219,489	43,845	4,690	1.8%	72.8%	82.6%	82.7%	72.7%	6929.8	±	3495.7	123.4	±	115.5	21.4	±	1.4	12.4	±	1.6
2011	33	252,749	232,572	20,177	3,741	1.5%	68.6%	89.6%	89.6%	68.6%	7659.1	±	3784.5	113.4	±	135.7	2.0	±	1.6	10.7	±	1.6
2012	33	262,296	204,833	57,463	5,373	2.0%	65.3%	88.4%	88.5%	65.2%	7948.4	±	4099.4	162.8	±	165.0	19.2	±	1.5	10.8	±	2.2
2013	31	225,389	169,323	56,066	5,041	2.2%	66.0%	89.3%	89.5%	65.8%	7270.6	±	3668.7	162.6	±	184.9	7.1	±	1.9	3.5	±	1.8
2014	31	218,341	207,921	10,420	5,190	2.4%	78.8%	93.2%	93.5%	78.5%	7043.3	±	3028.2	167.4	±	208.3	16.1	±	1.7	22.3	±	2.1
2015	31	247,419	211,110	36,309	5,146	2.1%	86.2%	89.6%	89.6%	86.2%	7981.3	±	3681.5	166.0	±	120.2	18.6	±	1.1	21.0	±	2.5
2016	29	237,275	204,648	32,627	4,800	2.0%	62.9%	81.6%	81.8%	62.7%	8181.9	±	4204.7	165.5	±	148.4	22.3	±	0.9	11.3	±	1.3
2017	14	100,999	100,140	859	1,485	1.5%	54.1%	87.0%	87.2%	53.9%	7214.2	±	3921.2	106.1	±	117.5	0.2	±	1.2	1.4	±	1.3
2018	19	157,551	143,196	14,355	3,206	2.0%	79.5%	82.4%	82.4%	79.5%	8292.2	±	3523.5	168.7	±	160.4	5.2	±	1.7	0.2	±	2.2
2019	14	105,547	97,185	8,362	1,319	1.2%	83.5%	86.7%	86.7%	83.5%	7539.1	±	3919.8	94.2	±	88.4	0.7	±	1.2	4.6	±	2.0
2020	13	110,426	70,789	39,637	1,692	1.5%	68.3%	72.2%	72.2%	68.3%	8494.3	±	3642.3	130.2	±	118.2	19.9	±	1.5	3.4	±	2.2

Notes

• lmk!

4.7.2 Figure 1

Let’s first download the map

Code

# with or without bathymetry
with_bathy <- F

# if bathy
if (with_bathy) {
  # check if background_ggoceanmaps exist
  if (file.exists("../export/background_ggoceanmap_zoom_out.rds")) {
    trip_zoom_out <- readRDS("../export/background_ggoceanmap_zoom_out.rds")
  } else {
    # using ggOceanMaps
    trip_zoom_out <- basemap(
      limits = c(170, -110, 30, 59),
      # limits = c(-155, -110, 30, 60),
      bathymetry = TRUE,
      shapefiles = "Arctic",
      rotate = TRUE,
      grid.col = NA
    )

    # Make the graticules:
    lims <- attributes(trip_zoom_out)$limits
    graticule <- sf::st_graticule(
      c(lims[1], lims[3], lims[2], lims[4]),
      crs = attributes(trip_zoom_out)$proj,
      lon = seq(-180, 180, 45),
      lat = seq(-90, 90, 10)
    )

    # Plot
    trip_zoom_out <- trip_zoom_out +
      geom_sf(data = graticule, color = "grey50")

    # reorder layers
    trip_zoom_out$layers <- trip_zoom_out$layers[c(1, 3, 2)]

    # save result
    saveRDS(trip_zoom_out, "../export/background_ggoceanmap_zoom_out.rds")
  }
} else {
  # using ggOceanMaps
  trip_zoom_out <- basemap(
    limits = c(170, -110, 20, 59),
    # limits = c(-155, -110, 30, 60),
    # land.col = 'black',
    # land.col = '#cdaa6d',
    bathymetry = FALSE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip_zoom_out)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip_zoom_out)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip_zoom_out <- trip_zoom_out +
    geom_sf(data = graticule, color = "grey70", linewidth = 0.3) +
    theme(
      panel.background = element_rect(fill = "#4292c6"),
      panel.ontop = FALSE
    )

  # reorder layers
  trip_zoom_out$layers <- trip_zoom_out$layers[c(2, 1)]
}

Tip

This part is mostly based on https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html. The author of the package confirms that his kernel density estimation is calculated using “meaningful” units such as UTM ¹, since the input is in a simple feature format.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat)) %>%
  # select only the required columns
  select(lat = Lat, long = Lon, id, DiveTypeName) %>%
  # create col to nicely display dive type
  mutate(origin = paste(DiveTypeName, "dives"))

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# make it's group_by origin
df_kernel_dens_sf <- group_by(df_kernel_dens_sf, origin)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(
  df_kernel_dens_sf,
  # H = diag(c(
  #   MASS::bandwidth.nrd(
  #     sf::st_coordinates(df_kernel_dens_sf)[, 1]
  #   ),
  #   MASS::bandwidth.nrd(
  #     sf::st_coordinates(df_kernel_dens_sf)[, 2]
  #   )
  # ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
to_plot <- trip_zoom_out +
  # rename axis
  labs(x = "Longitude", y = "Latitude") +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probability contours") +
  # no display alpha
  guides(alpha = "none") +
  # facet by origin
  facet_wrap(. ~ origin)

# reorder layers
to_plot$layers <- if (with_bathy) {
  to_plot$layers[c(1, 2, 4, 3)]
} else {
  to_plot$layers[c(1, 3, 2)]
}

# display
to_plot

Figure 7: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probability contours represent different levels of probability density.

Notes

• lmk!

4.7.3 Figure 2

Code

# the last plot
prop_at_sea_trip <- data_dive %>%
  group_by(id, trip) %>%
  # time difference between the first date and the others (telling R to take diff b/w first date and all other dates)
  mutate(nb_days_departure = trunc(as.numeric(difftime(
    date, first(date),
    units = "days"
  )))) %>%
  # group by day/dive_type/seal
  group_by(nb_days_departure, DiveTypeName, id, trip) %>%
  # count number of dives
  summarise(nb_daily_dives_type = n(), .groups = "drop") %>%
  # group by day/seal
  group_by(nb_days_departure, id) %>%
  # count the total number of dives per day/seal
  mutate(nb_daily_dives = sum(nb_daily_dives_type)) %>%
  # order by seals/date
  arrange(id, nb_days_departure, trip) %>%
  # perform our calculation per seal
  group_by(id, trip) %>%
  # calculate of the Proportion of time since departure
  mutate(perc_time_at_sea = round(
    nb_days_departure / max(nb_days_departure) * 100, 1
  )) %>%
  # filter on benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # add class of Proportion of time at sea
  mutate(day_class = cut(perc_time_at_sea, seq(0, 100, 5),
    include.lowest = T
  )) %>%
  # by class of Proportion of time at sea
  group_by(day_class, trip) %>%
  # calculate the proportion of daily benthic dives
  summarise(
    nb_benthic_class = sum(nb_daily_dives_type),
    nb_total_class = sum(nb_daily_dives),
    prop_class = nb_benthic_class / nb_total_class
  ) %>%
  # plot
  ggplot(
    .,
    aes(x = day_class, y = prop_class)
  ) +
  geom_bar(stat = "identity", fill = "#008c49") +
  guides(x = guide_axis(angle = 45)) +
  scale_y_continuous(
    limits = c(0, 1),
    # breaks=seq(0,0.6,0.1),
    labels = function(x) {
      scales::percent(abs(x), 1)
    }
  ) +
  facet_wrap(. ~ trip) +
  # annotate("text", x = 5.3, y = 0.4, label = "Departure") +
  # annotate("text", x = 15, y = 0.5, label = "Arrival") +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(4),
  #     x2 = c(1),
  #     y1 = c(0.4),
  #     y2 = c(0.33)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = 0.3
  # ) +
  # geom_curve(
  #   data =   tibble(
  #     x1 = c(16),
  #     x2 = c(19.7),
  #     y1 = c(0.5),
  #     y2 = c(0.57)
  #   ),
  #   aes(
  #     x = x1,
  #     y = y1,
  #     xend = x2,
  #     yend = y2
  #   ),
  #   arrow = arrow(length = unit(0.08, "inch")),
  #   linewidth = 0.5,
  #   color = "gray20",
  #   curvature = -0.56
  # ) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(arrow = arrow(
      angle = 20,
      length = unit(.10, "inches"),
      type = "closed"
    ))
  ) +
  labs(
    x = "Time spent at sea (%)",
    y = "Proportion of benthic dives"
  )

Code

# display
prop_at_sea_trip

Figure 8: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea during post breeding and post molting trip. Time spent at sea is classified into intervals of five percent for graphical representation.

4.8 Extra

Code

# import the csv export
data_plot <- readr::read_csv("../data/benthic_dive_extract.csv")

# plot
data_plot %>%
  ggplot(aes(x = time, y = CorrectedDepth)) +
  geom_path() +
  scale_y_reverse() +
  theme_minimal() +
  theme(
    axis.line.y = element_line(arrow = arrow(
      type = "closed",
      ends = "first",
      length = unit(0.1, "inches")
    )),
    axis.line.x = element_line(arrow = arrow(
      type = "closed",
      ends = "last",
      length = unit(0.1, "inches")
    ))
  ) +
  labs(
    x = "Time",
    y = "Depth (m)"
  ) +
  # geom_text(data = data.frame(
  #          xpos = median(data_plot$time),
  #          ypos = 30,
  #          hjustvar = 0,
  #          vjustvar = 0,
  #          annotateText = paste0("Corner index: ", round(unique(data_plot$corner_index),2),
  #                       "\n",
  #                       "Benthix index: ", round(unique(data_plot$benthic_index)),2)),
  #           aes(x=xpos,y=ypos,hjust=hjustvar,
  #               vjust=vjustvar,label=annotateText)) +
  scale_x_datetime(position = "top") +
  labs(subtitle = "ID: 2011034 - Dive: 29")

Footnotes

1. • https://gis.stackexchange.com/questions/64638/how-to-fill-in-the-parameters-for-kernel-density-estimation
   • https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1020&context=siv
   • https://gis.stackexchange.com/questions/50692/choosing-projection-for-rasterization-of-random-latitude-and-longitude-data-in-n/50791#50791